# ------------------------------------------------------------------------------#
#	Paper:			Local Effects of Foreign Direct Investment
#	Authors:		Tobias Rommel, Tabea Palmtag, Luca Messerschmidt
#	Version:		January 29, 2025
#	Data:			  FDI projects, nightlights data, V-Dem data
#	Task:			  Treatment matching analysis
# ------------------------------------------------------------------------------#


# Matching Methods for Causal Inference with Time-Series Cross-Sectional Data (Kosuke Imai, In Song Kim, and Erik H. Wang)
# Please note that this code was running for several hours on a cloud computer of the Leibniz Supercomputing Centre. 
# Errors in estimating these models with a normal computer might occur and the process might be aborted due to memory exhaustion.


## Load Packages
library(sf)
library(tidyverse)
library(countrycode)
library(dplyr)
library(stars)
library(haven)
library(readr)
library(PanelMatch)


## Directory
rm(list=ls())
setwd("/Users/tobiasrommel/Downloads/RPM_FDIgrowth_Replication/")


# ------------------------------------------------------------------------------#


### Treatment Matching Graphs, based on previous estimates (see below)


## Figure 5
dat <- readRDS("FDIgrowth_both_all.Rds") %>% slice(1:3) %>% dplyr::mutate(year = row_number()-1, `Grid Size`= "10km") %>% rename("low"=3,"up"=4) 
# Plot
ggplot(dat, aes(x = year, y = estimate, color = `Grid Size`)) +  
  geom_errorbar(
    aes(ymin = low, ymax = up),  
    position = position_dodge(0.3), width = 0.2) +
  geom_point(position = position_dodge(0.3)) +  
  scale_x_continuous(breaks = c(0, 1, 2)) +  
  ylab("Coefficient Estimate") +
  xlab("") +
  ggtitle("Nightlight Intensity, Unconditional (10x10km)") +  
  theme_light() +
  geom_hline(yintercept = 0, color = "darkgrey", size = 0.5) +
  theme(legend.position = "none")
ggsave("Figure5.eps", width = 6, height = 4, dpi = 300)


## Figure 6
dat <- readRDS("FDIspatial_both_all.Rds") %>% slice(1:3) %>% dplyr::mutate(year = row_number()-1, `Grid Size`= "25km") %>% rename("low"=3,"up"=4)
# Plot
ggplot(dat, aes(x = year, y = estimate, color = `Grid Size`)) +  
  geom_errorbar(
    aes(ymin = low, ymax = up),  
    position = position_dodge(0.3), width = 0.2) +
  geom_point(position = position_dodge(0.3)) + 
  scale_x_continuous(breaks = c(0, 1, 2)) +  
  ylab("Coefficient Estimate") +
  xlab("") +
  ggtitle("Nightlight Inequality, Unconditional (25x25km)") +  
  theme_light() +
  geom_hline(yintercept = 0, color = "darkgrey", size = 0.5) +
  theme(legend.position = "none")
ggsave("Figure6.eps", width = 6, height = 4, dpi = 300)


## Figure A-3
aut <- readRDS("FDIgrowth_aut_all.Rds") %>% slice(1:3) %>% dplyr::mutate(year = row_number()-1,  Type = "Autocracies") %>% rename("low"=3,"up"=4)
dem <- readRDS("FDIgrowth_dem_all.Rds") %>% slice(1:3) %>% dplyr::mutate(year = row_number()-1,  Type = "Democracies") %>% rename("low"=3,"up"=4)
dat <- bind_rows(aut, dem)

# Plot
ggplot(dat, aes(x = year, y = estimate, color = Type)) +  # Use Type for different colors
  geom_errorbar(
    aes(ymin = low, ymax = up),
    position = position_dodge(0.3), width = 0.2
  ) +
  geom_point(position = position_dodge(0.3)) +
  scale_x_continuous(breaks = c(0, 1, 2)) +
  ylab("Coefficient Estimate") +
  xlab("") +
  ggtitle("Nightlight Intensity, by Regime Type (10x10km)") +
  theme_light() +
  geom_hline(yintercept = 0, color = "darkgrey", size = 0.5) +
  theme(legend.position = "bottom", legend.title = element_blank()) 
ggsave("FigureA3.eps", width = 6, height = 4, dpi = 300)


## Figure A-4
aut <- readRDS("FDIspatial_aut_all.Rds") %>% slice(1:3) %>% dplyr::mutate(year = row_number()-1,  Type = "Autocracies") %>% rename("low"=3,"up"=4)
dem <- readRDS("FDIspatial_dem_all.Rds") %>% slice(1:3) %>% dplyr::mutate(year = row_number()-1,  Type = "Democracies") %>% rename("low"=3,"up"=4)
dat <- bind_rows(aut, dem)
# Plot
ggplot(dat, aes(x = year, y = estimate, color = Type)) +  # Use Type for different colors
  geom_errorbar(
    aes(ymin = low, ymax = up),
    position = position_dodge(0.3), width = 0.2
  ) +
  geom_point(position = position_dodge(0.3)) +
  scale_x_continuous(breaks = c(0, 1, 2)) +
  ylab("Coefficient Estimate") +
  xlab("") +
  ggtitle("Nightlight Inequality, by Regime Type (25x25km)") +
  theme_light() +
  geom_hline(yintercept = 0, color = "darkgrey", size = 0.5) +
  theme(legend.position = "bottom", legend.title = element_blank()) 
ggsave("FigureA4.eps", width = 6, height = 4, dpi = 300)



# ------------------------------------------------------------------------------#


### Treatment Matching Estimations (will take several hours, might run out of memory!)


## Load 25x25km Data (adapt for 10x10km analysis)
dt_matching_all <- read_dta("RPM_FDIgrowth_Grid25km.dta")
saveRDS(dt_matching_all,("dt_matching_all.Rds"))

## Nightlight Intensity

# Overall
dt_matching<- readRDS(paste0("dt_matching_all.Rds")) 
dt_matching$fdi_treat<- as.integer(dt_matching$fdi_treat)
dt_matching$grid<- as.integer(dt_matching$grid)
dt_matching$lights_mean<- as.integer(dt_matching$lights_mean)
dt_matching$hyde_mean<- as.integer(dt_matching$hyde_mean)
dt_matching$year<-as.integer(dt_matching$year)
dt_matching$cname<-as.factor(dt_matching$cname)
dt_matching <- dt_matching[!is.na(dt_matching$grid),]
dt_matching<-as.data.frame(dt_matching)
dt_matching<-dt_matching%>%dplyr::rename("id"="grid") %>% dplyr::select(1:22)
PM.match <- PanelMatch(lag = 3, time.id = "year", unit.id = "id",
                       treatment = "fdi_treat", refinement.method = "ps.match",
                       data = dt_matching, match.missing = TRUE,
                       covs.formula = ~ cname + I(lag(hyde_mean, 1:2))+ I(lag(lights_mean, 1:2)), 
                       size.match = 5, qoi = "att",
                       outcome.var = "lights_mean", lead = 0:10, forbid.treatment.reversal = FALSE)
PE.results <- PanelEstimate(sets = PM.match, data = dt_matching)
a<- summary(PE.results)
PM.results<- as.data.frame(a$summary)
saveRDS(PM.results, paste0("FDIgrowth_both_all.Rds"))

# Democracy
dt_matching<- readRDS(paste0("dt_matching_all.Rds")) %>%dplyr::filter(democracy==1)
dt_matching$fdi_treat<- as.integer(dt_matching$fdi_treat)
dt_matching$grid<- as.integer(dt_matching$grid)
dt_matching$lights_mean<- as.integer(dt_matching$lights_mean)
dt_matching$hyde_mean<- as.integer(dt_matching$hyde_mean)
dt_matching$year<-as.integer(dt_matching$year)
dt_matching$cname<-as.factor(dt_matching$cname)
dt_matching <- dt_matching[!is.na(dt_matching$grid),]
dt_matching<-as.data.frame(dt_matching)
dt_matching<-dt_matching%>%dplyr::rename("id"="grid") %>% dplyr::select(1:22)
PM.match <- PanelMatch(lag = 3, time.id = "year", unit.id = "id",
                       treatment = "fdi_treat", refinement.method = "ps.match",
                       data = dt_matching, match.missing = TRUE,
                       covs.formula = ~ cname + I(lag(hyde_mean, 1:2))+ I(lag(lights_mean, 1:2)), 
                       size.match = 5, qoi = "att",
                       outcome.var = "lights_mean", lead = 0:10, forbid.treatment.reversal = FALSE)
PE.results <- PanelEstimate(sets = PM.match, data = dt_matching)
a<- summary(PE.results)
PM.results<- as.data.frame(a$summary)
saveRDS(PM.results, paste0("FDIgrowth_dem_all.Rds"))

# Autocracy
dt_matching<- readRDS(paste0("dt_matching_all.Rds")) %>%dplyr::filter(democracy==0)
dt_matching$fdi_treat<- as.integer(dt_matching$fdi_treat)
dt_matching$grid<- as.integer(dt_matching$grid)
dt_matching$lights_mean<- as.integer(dt_matching$lights_mean)
dt_matching$hyde_mean<- as.integer(dt_matching$hyde_mean)
dt_matching$year<-as.integer(dt_matching$year)
dt_matching$cname<-as.factor(dt_matching$cname)
dt_matching <- dt_matching[!is.na(dt_matching$grid),]
dt_matching<-as.data.frame(dt_matching)
dt_matching<-dt_matching%>%dplyr::rename("id"="grid") %>% dplyr::select(1:22)
PM.match <- PanelMatch(lag = 3, time.id = "year", unit.id = "id",
                       treatment = "fdi_treat", refinement.method = "ps.match",
                       data = dt_matching, match.missing = TRUE,
                       covs.formula = ~ cname + I(lag(hyde_mean, 1:2))+ I(lag(lights_mean, 1:2)), 
                       size.match = 5, qoi = "att",
                       outcome.var = "lights_mean", lead = 0:10, forbid.treatment.reversal = FALSE)
PE.results <- PanelEstimate(sets = PM.match, data = dt_matching)
a<- summary(PE.results)
PM.results<- as.data.frame(a$summary)
saveRDS(PM.results, paste0("FDIgrowth_aut_all.Rds"))


## Nightlight Inequality

# Overall
dt_matching<- readRDS(paste0("dt_matching_all.Rds")) 
dt_matching$fdi_treat<- as.integer(dt_matching$fdi_treat)
dt_matching$grid<- as.integer(dt_matching$grid)
dt_matching$lights_diff<- as.integer(dt_matching$lights_diff)
dt_matching$hyde_mean<- as.integer(dt_matching$hyde_mean)
dt_matching$year<-as.integer(dt_matching$year)
dt_matching$cname<-as.factor(dt_matching$cname)
dt_matching <- dt_matching[!is.na(dt_matching$grid),]
dt_matching<-as.data.frame(dt_matching)
dt_matching<-dt_matching%>%dplyr::rename("id"="grid") %>% dplyr::select(1:22)
PM.match <- PanelMatch(lag = 3, time.id = "year", unit.id = "id",
                       treatment = "fdi_treat", refinement.method = "ps.match",
                       data = dt_matching, match.missing = TRUE,
                       covs.formula = ~ cname + I(lag(hyde_mean, 1:2))+ I(lag(lights_diff, 1:2)), 
                       size.match = 5, qoi = "att",
                       outcome.var = "lights_diff", lead = 0:10, forbid.treatment.reversal = FALSE)
PE.results <- PanelEstimate(sets = PM.match, data = dt_matching)
a<- summary(PE.results)
PM.results<- as.data.frame(a$summary)
saveRDS(PM.results, paste0("FDIspatial_both_all.Rds"))

# Democracy
dt_matching<- readRDS(paste0("dt_matching_all.Rds")) %>%dplyr::filter(democracy==1) 
dt_matching$fdi_treat<- as.integer(dt_matching$fdi_treat)
dt_matching$grid<- as.integer(dt_matching$grid)
dt_matching$lights_diff<- as.integer(dt_matching$lights_diff)
dt_matching$hyde_mean<- as.integer(dt_matching$hyde_mean)
dt_matching$year<-as.integer(dt_matching$year)
dt_matching$cname<-as.factor(dt_matching$cname)
dt_matching <- dt_matching[!is.na(dt_matching$grid),]
dt_matching<-as.data.frame(dt_matching)
dt_matching<-dt_matching%>%dplyr::rename("id"="grid") %>% dplyr::select(1:22)
PM.match <- PanelMatch(lag = 3, time.id = "year", unit.id = "id",
                       treatment = "fdi_treat", refinement.method = "ps.match",
                       data = dt_matching, match.missing = TRUE,
                       covs.formula = ~ cname + I(lag(hyde_mean, 1:2))+ I(lag(lights_diff, 1:2)), 
                       size.match = 5, qoi = "att",
                       outcome.var = "lights_diff", lead = 0:10, forbid.treatment.reversal = FALSE)
PE.results <- PanelEstimate(sets = PM.match, data = dt_matching)
a<- summary(PE.results)
PM.results<- as.data.frame(a$summary)
saveRDS(PM.results, paste0("FDIspatial_dem_all.Rds"))

# Autocracy
dt_matching<- readRDS(paste0("dt_matching_all.Rds")) %>%dplyr::filter(democracy==0)
dt_matching$fdi_treat<- as.integer(dt_matching$fdi_treat)
dt_matching$grid<- as.integer(dt_matching$grid)
dt_matching$lights_diff<- as.integer(dt_matching$lights_diff)
dt_matching$hyde_mean<- as.integer(dt_matching$hyde_mean)
dt_matching$year<-as.integer(dt_matching$year)
dt_matching$cname<-as.factor(dt_matching$cname)
dt_matching <- dt_matching[!is.na(dt_matching$grid),]
dt_matching<-as.data.frame(dt_matching)
dt_matching<-dt_matching%>%dplyr::rename("id"="grid") %>% dplyr::select(1:22)
PM.match <- PanelMatch(lag = 3, time.id = "year", unit.id = "id",
                       treatment = "fdi_treat", refinement.method = "ps.match",
                       data = dt_matching, match.missing = TRUE,
                       covs.formula = ~ cname + I(lag(hyde_mean, 1:2))+ I(lag(lights_diff, 1:2)), 
                       size.match = 5, qoi = "att",
                       outcome.var = "lights_diff", lead = 0:10, forbid.treatment.reversal = FALSE)
PE.results <- PanelEstimate(sets = PM.match, data = dt_matching)
a<- summary(PE.results)
PM.results<- as.data.frame(a$summary)
saveRDS(PM.results, paste0("FDIspatial_aut_all.Rds"))





# ------------------------------------------------------------------------------#
